#include "som.h"

//-----------------------------------------------------------------------------------------------------------
SOM::SOM(){
	//printf("Entra al constructor. \n");
	loadConfig();
	init();
}

SOM::SOM(vector< vector<double> >& v){
	loadConfig();
	init();
	this->inputs.resize(v.size());
	for (unsigned int i=0;i<v.size();i++){
		this->inputs[i].inp.resize(v[i].size());
		for (unsigned int j=0;j<v[i].size();j++){
			this->inputs[i].inp[j] = v[i][j];
		}
	}	
	//normalize();
}

void SOM::init(){
	//Generamos la malla (vNeurons), con el mapeo nodo-conectividades (vIcone) y pesos auxiliares.
	int i,j,k,numb;
	this->vIcone.resize(this->x*this->y);
	this->vNeurons.resize(this->x*this->y);
	this->vPos.resize(this->x*this->y);
	srand(time(NULL)); //semilla para el random
	if (this->x == 1){ //Si x==1
		for (j=0;j<this->y;j++){
			numb = j;
			this->vPos[numb].x = 1;
			this->vPos[numb].y = j;
			if (!j){
				this->vIcone[numb].conec.resize(1);
				this->vIcone[numb].conec[0] = numb+1;
			}
			else{
				if (j == (this->y-1)){
					this->vIcone[numb].conec.resize(1);
					this->vIcone[numb].conec[0] = numb-1;
				}
				else{
					this->vIcone[numb].conec.resize(2);
					this->vIcone[numb].conec[0] = numb-1;
					this->vIcone[numb].conec[1] = numb+1;
				}
			}/*
			cout<<"Neurona "<<j<<": ";
			for (int k=0;k<vIcone[numb].conec.size();k++){
				cout<<" "<<this->vIcone[numb].conec[k];
			}
			cout<<endl;*/
			vNeurons[numb].weights.resize(this->dInit); //Cargo pesos iniciales
			for (k=0;k<this->dInit;k++) vNeurons[numb].weights[k] = randomWeight();
		}
	}
	else{
		if (this->y == 1){ //Si y==1
			for (i=0;i<this->x;i++){
				numb = i;
				this->vPos[numb].x = i;
				this->vPos[numb].y = 1;
				if (!i){
					this->vIcone[numb].conec.resize(1);
					this->vIcone[numb].conec[0] = numb+1;
				}
				else{
					if (i == (this->x-1)){
						this->vIcone[numb].conec.resize(1);
						this->vIcone[numb].conec[0] = numb-1;
					}
					else{
						this->vIcone[numb].conec.resize(2);
						this->vIcone[numb].conec[0] = numb-1;
						this->vIcone[numb].conec[1] = numb+1;
					}
				}/*
				cout<<"Neurona "<<i<<": ";
				for (int k=0;k<vIcone[numb].conec.size();k++){
					cout<<" "<<this->vIcone[numb].conec[k];
				}
				cout<<endl;*/
				vNeurons[numb].weights.resize(this->dInit); //Cargo pesos iniciales
				for (k=0;k<this->dInit;k++) vNeurons[numb].weights[k] = randomWeight();
			}
		}
		else{
			for (i=0;i<this->x;i++){
				for (j=0;j<this->y;j++){
					numb = i*this->y+j;
					this->vPos[numb].x = i;
					this->vPos[numb].y = j;
					if (!i){//Lateral izquierda
						if (!j){
							//Esquina inferior izquierda
							this->vIcone[numb].conec.resize(2);
							this->vIcone[numb].conec[0] = numb+1;
							this->vIcone[numb].conec[1] = numb+this->y;
						}
						else{ 
							if (j<(y-1)){
								//medio
								this->vIcone[numb].conec.resize(3);
								this->vIcone[numb].conec[0] = numb-1;
								this->vIcone[numb].conec[1] = numb+1;
								this->vIcone[numb].conec[2] = numb+this->y;
							}
							else{
								//arriba
									this->vIcone[numb].conec.resize(2);
									this->vIcone[numb].conec[0] = numb-1;
									this->vIcone[numb].conec[1] = numb+this->y;
							}
						}  
					}//i!=0
					else{ 
						if (!j){//Inferior 1,0,2.0,...
							if (i<(this->x-1)){
								this->vIcone[numb].conec.resize(3);
								this->vIcone[numb].conec[0] = numb-this->y;
								this->vIcone[numb].conec[1] = numb+1;
								this->vIcone[numb].conec[2] = numb+this->y;
							}
							else{//Esquina inferior derecha
								this->vIcone[numb].conec.resize(2);
								this->vIcone[numb].conec[0] = numb+1;
								this->vIcone[numb].conec[1] = numb-this->y;
							}
						}//if!j
						else{//j!=0
							//Estoy dentro de la grilla incluyendo lateral izquierda y parte superior
							if ((i<(this->x-1)) && (j<(this->y-1))){
								//Estoy dentro de la grilla sin incluir lateral izquierda y parte superior
								this->vIcone[numb].conec.resize(4);
								this->vIcone[numb].conec[0] = numb-1;
								this->vIcone[numb].conec[1] = numb-this->y;
								this->vIcone[numb].conec[2] = numb+1;
								this->vIcone[numb].conec[3] = numb+this->y;
							}
							else{
								if ((i==(this->x-1)) && (j==(this->y-1))){
									//Esquina superior derecha (ultimo nodo en la numeracion)
									this->vIcone[numb].conec.resize(2);
									this->vIcone[numb].conec[0] = numb-this->y;
									this->vIcone[numb].conec[1] = numb-1;
								}
								else{
									if ((i==(this->x-1)) && (j<(this->y-1))){
										//Lateral derecha
										this->vIcone[numb].conec.resize(3);
										this->vIcone[numb].conec[0] = numb-1;
										this->vIcone[numb].conec[1] = numb-this->y;
										this->vIcone[numb].conec[2] = numb+1;
									}
									else{
										//Parte superior
										if ((i<(this->x-1)) && (j==(this->y-1))){
										this->vIcone[numb].conec.resize(3);
										this->vIcone[numb].conec[0] = numb-this->y;
										this->vIcone[numb].conec[1] = numb-1;
										this->vIcone[numb].conec[2] = numb+this->y;
									}						
								}
							}
						}
					}	
				}	/*
				cout<<"Nodo número: "<<numb<<" y sus conectividades: ";
				for	(k=0;k<this->vIcone[numb].conec.size();k++){
					cout<<this->vIcone[numb].conec[k]<<" ";
				}
				cout<<endl;*/
				vNeurons[numb].weights.resize(this->dInit); //Cargo pesos iniciales
				for (k=0;k<this->dInit;k++) vNeurons[numb].weights[k] = randomWeight();
				}//for j
			}//for i
		}
	}
}
//-----------------------------------------------------------------------------------------------------------
/*Proceso de competición, fórmula pag 470
		i(x)=argmin(x-wj)
*/
int SOM::competency(int& k){
	unsigned int i,j;
	double min = 100000,ix,winner;
	k = rand() % this->inputs.size(); //Seleccionamos una entrada aleatoria -> mejorarlo con algo que no permita que en la misma epoca tome dos iguales!
	for (i=0;i<this->vNeurons.size();i++){//Para cada neurona
		ix = 0;
		//para cada dimension (en la entrada k-esima)
		for (j=0;j<this->dInit;j++){ 
			ix += pow(this->inputs[k].inp[j]-this->vNeurons[i].weights[j],2);
		}
		ix = sqrt(ix);
		if (ix<min){
			min = ix;
			winner = i; //Winning neuron
		}
	}
	//cout<<winner<<endl;	cin.get();
	return winner;
}
//-----------------------------------------------------------------------------------------------------------
void SOM::loadConfig(){
	//Cargamos las neuronas y sus conectividades
	ifstream archi;
	archi.open("config.txt");
	if(archi.is_open()){
		string str_line;
		char* line;

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->epOrder = atoi(line);	

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->epTrans = atoi(line);	

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->epConv = atoi(line);	

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->eta = atof(line);	
		//printf("%f \n",this->eta0);

		this->eta0 = this->eta;

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->dInit = atof(line);
		//printf("%f \n",this->dInit);

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->numbInputs = atoi(line);
		//printf("%i \n",this->numbInputs);

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->x = atoi(line);
		//printf("%i \n",this->x);

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->y = atoi(line);
		//printf("%i \n",this->y);

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->style = atoi(line);

		getline(archi,str_line);
		getline(archi,str_line);
		line = &str_line[0];
		this->stPrint = atoi(line);

		//this->tau2 = epocas;
		//printf("%i \n",this->tau2);

		//this->sigma0 = (this->x > this->y) ? (double)this->x/2.0 : (double)this->y/2.0;

		this->sigma = (this->x > this->y) ? (double)this->x/2.0 : (double)this->y/2.0; //Ag

		this->sigma0 = this->sigma;

		//this->tau1 = (double)(this->epocas)/log(this->sigma0);
		//this->tau1 = 8*this->tau1;
		//this->tau1 = 100;
		//printf("%f \n",this->tau1);

		//this->eta = 0.1;	

		archi.close();
	}
	else cout<<"Imposible cargar archivo de configuraciones."<<endl;
}
//-----------------------------------------------------------------------------------------------------------
float uniform(){
	return ((double)rand() / ((double)(RAND_MAX)+(double)(1)));
}

double SOM::randomWeight(){
	double f = uniform();
	return (double)f - 0.5;
}
//-----------------------------------------------------------------------------------------------------------
/*Actualizamos los pesos de la ganadora y su vecindario. Fórmula 9.13 pág. 473*/
void SOM::Neighborhood(int& winner,int& n,int& k,int phase){
	double dist,Hji;
	int tam; 
	switch (phase){
		case 0 : { //Ordenamiento
			for (unsigned int i=0;i<this->vNeurons.size();i++){
				dist = distance(this->vPos[i].x,this->vPos[i].y,this->vPos[winner].x,this->vPos[winner].y); //Distancia
				Hji = exp(-dist/(2.0*(this->sigma*this->sigma)));
				//cout<<"Hji : "<<Hji<<endl;
				for (unsigned int j=0;j<this->vNeurons[i].weights.size();j++){
					this->vNeurons[i].weights[j] = this->vNeurons[i].weights[j]+this->eta*Hji*(this->inputs[k].inp[j]-this->vNeurons[i].weights[j]);
					//cout<<this->vNeurons[i].weights[j]<<endl;
				}
			}
			//cin.get();
			break;
		}
		case 1 : { //Transicion
			tam = this->epTrans*this->inputs.size();
			this->eta = ((0.1-this->eta0)/tam)*n+this->eta0;
			//this->eta = this->eta0*exp(-n/tam);
			//cout<<this->eta<<endl;	
			this->sigma = ((1 - this->sigma0)/tam)*n+this->sigma0;
			//this->sigma = this->sigma0*exp(-n/tam);
			for (unsigned int i=0;i<this->vNeurons.size();i++){
				//Observar como eta y sigma decrecen linealmente	
				dist = distance(this->vPos[i].x,this->vPos[i].y,this->vPos[winner].x,this->vPos[winner].y); //Distancia
				Hji = exp(-dist/(2.0*(this->sigma*this->sigma))); //como sigma decrece, el argumento de la exp es mas chico ->Hji decrece
				//cout<<Hji<<endl; cin.get();
				for (unsigned int j=0;j<this->vNeurons[i].weights.size();j++){
					this->vNeurons[i].weights[j] = this->vNeurons[i].weights[j]+this->eta*Hji*(this->inputs[k].inp[j]-this->vNeurons[i].weights[j]);
					//cout<<this->vNeurons[i].weights[j]<<endl;
				}
			}			
			break;
		}
		case 2 : { //Convergencia
				tam = this->epConv*this->inputs.size();
				this->eta = ((0.01-0.1)/tam)*n+0.1;
				//cout<<this->eta<<endl;
				for (unsigned int j=0;j<this->vNeurons[winner].weights.size();j++){
						this->vNeurons[winner].weights[j] = this->vNeurons[winner].weights[j]+this->eta*(this->inputs[k].inp[j]-this->vNeurons[winner].weights[j]);
				}
				//cout<<this->vNeurons[i].weights[j]<<endl;
			break;
		}
	}
}
//-----------------------------------------------------------------------------------------------------------
/*Es la distancia que esta en el libro pag 472, i es la ganadora*/
double SOM::distance(int& rjx,int& rjy,int& rix,int& riy){
	double dist = (rjx-rix)*(rjx-rix)+(rjy-riy)*(rjy-riy);	
	//cout<<dist<<endl;
	return dist;
	//return abs(rjx-rix)+abs(rjy-riy);
}
//-----------------------------------------------------------------------------------------------------------
//Global para poner replotear en la misma instancia!
GNUplot gp; 

/*Para printear el mapa*/
void SOM::print(bool classif){
	int i,len = this->vNeurons.size();
	//Borro las referencias de GNUPlot	
	stringstream ss3;
	ss3.flags(std::ios::fixed);
	ss3<<"set nokey";
	gp(ss3.str());
	stringstream ss2;
	ss2.flags(std::ios::fixed);
	if(!classif){
		//Guardo los pesos de cada Neurona por separado
		for (i=0;i<len;i++){
			stringstream ss1;
			ss1<<"neurWeights";
			ss1<<i<<".txt";
			saveWeightsFile(i,ss1.str());
		}
		ss2<<"plot [-1:1][-1:1] ";
		switch(style){
				case 0 : { 
					for (int i=0;i<(len-1);i++) ss2<<"'neurWeights"<<i<<".txt' with lines lt 1, ";
					ss2<<"'neurWeights"<<(len-1)<<".txt' with lines lt 1;";
					break;
				}
				case 1 : { 
					for (int i=0;i<(len-1);i++) ss2<<"'neurWeights"<<i<<".txt' with points lt 1, ";
					ss2<<"'neurWeights"<<(len-1)<<".txt' with points lt 1;";
					break;
				}
				case 2 : { 
					for (int i=0;i<(len-1);i++) ss2<<"'neurWeights"<<i<<".txt' with linespoints lt 1, ";
					ss2<<"'neurWeights"<<(len-1)<<".txt' with linespoints lt 1;";
					break;
				}
		}
	}
	else{
		vector<neurons> Class1,Class2;
		for (i=0;i<len;i++){
			if (!this->vClass[i]) Class1.push_back(vNeurons[i]);
			else Class2.push_back(vNeurons[i]);
		}		
		saveWeightsFile(Class1,"class1.txt");
		saveWeightsFile(Class2,"class2.txt");
		ss2<<"plot [-3:3][-3:3] 'class1.txt' with points lt 1, 'class2.txt' with points lt 3;";
	}	
	gp(ss2.str());
}
//-----------------------------------------------------------------------------------------------------------
void SOM::ordering_phase(bool classif){
	int winner,k;
	printf("Comienza la fase de ordenamiento. \n");
	for (int i=0;i<this->epOrder*this->inputs.size();i++){
		winner = competency(k); //Las neuronas compiten por una entrada
		Neighborhood(winner,i,k,0);
		if (!classif){
			if (i%(int)((this->epOrder*this->inputs.size())*0.03)==0) 
				print(0);//Dibujo
		}
	}
	//cout<<"Salio de ordenamiento"<<endl;
}
//-----------------------------------------------------------------------------------------------------------
void SOM::transition_phase(bool classif){
	int winner,k;
	printf("Comienza la fase de transición. \n");
	for (int i=0;i<this->epTrans*this->inputs.size();i++){
		winner = competency(k); //Las neuronas compiten por una entrada
		Neighborhood(winner,i,k,1);
		if (!classif){	  
			if (i%(int)((this->epTrans*this->inputs.size())*0.03)==0) 
				print(0);//Dibujo
		}
	}
	//cout<<"Salio de ordenamiento"<<endl;
}
//-----------------------------------------------------------------------------------------------------------
void SOM::convergence_phase(bool classif){
	int winner,k;
	printf("Comienza la fase de convergencia. \n");
	for (int i=0;i<this->epConv*this->inputs.size();i++){
		winner = competency(k); //Las neuronas compiten por una entrada
		Neighborhood(winner,i,k,2);
		if (!classif){		
			if (i%(int)((epConv*this->inputs.size())*0.03)==0) 
				print(0);//Dibujo
		}
	}
	//cout<<"Salio de convergencia"<<endl;
}
//-----------------------------------------------------------------------------------------------------------

#define PI 3,14159265

void SOM::generate_patterns(){
	//n==0 square patternes [-1:1][-1:1]
	//n==1 unit circle r=1 c=0
	//n==2 T
	double aux,r;
	this->inputs.resize(this->numbInputs);
	switch (this->stPrint){
		case 0 : {
				cout<<"Entradas en cuadrado unitario. \n";
				for (int i=0;i<this->numbInputs;i++){
					this->inputs[i].inp.resize(this->dInit);
					for (int j=0;j<this->dInit;j++){
						this->inputs[i].inp[j] = 2*randomWeight();
						//cout<<this->inputs[i].inp[j]<<endl;
					}
				}
				break;
		}
		case 1 : {
				cout<<"Entradas en circulo unitario. \n";
				for (int i=0;i<this->numbInputs;i++){
					this->inputs[i].inp.resize(this->dInit);
					aux = rand() % 360;
					r = (double)rand() / (double)RAND_MAX;
					this->inputs[i].inp[0] = r*cos((PI/180.0)*aux);
					this->inputs[i].inp[1] = r*sin((PI/180.0)*aux);
					//cout<<"Pasada ["<<i<<"] : x="<<this->inputs[i].inp[0]<<" y="<<this->inputs[i].inp[1]<<endl;
				}
				break;
		}
		case 2 : {
				cout<<"Entradas en T. \n";
				for (int i=0;i<this->numbInputs;i++){
					this->inputs[i].inp.resize(this->dInit);
					this->inputs[i].inp[0] = 2.0*randomWeight();
					if (fabs(this->inputs[i].inp[0])<=0.25) this->inputs[i].inp[1] = 2.0*randomWeight();
					else this->inputs[i].inp[1] = (randomWeight())/2.0 + 0.75;
				}
				break;
		}
		default : {
			printf("La opción %i no es válida.",this->stPrint);			
			break;
		}
	}
	//normalize(); //Para evitar errores numéricos
}
//-----------------------------------------------------------------------------------------------------------
/*Para ver si las conectividades son correctas*/
void SOM::saveConectivityFile(int& i,string file){
	ofstream archi(file.c_str(),ios::trunc);
	int xa,ya,xb,yb;
	if(archi.is_open()){	
		for (unsigned int n=0;n<vIcone[i].conec.size();n++){
				xa = vPos[i].x;
				ya = vPos[i].y;
				xb = vPos[vIcone[i].conec[n]].x;
				yb = vPos[vIcone[i].conec[n]].y;
				archi<<xa<<"       "<<ya<<endl;
				archi<<xb<<"       "<<yb<<endl;
		}
		archi.close();
	}
}
//-----------------------------------------------------------------------------------------------------------
/*Este si es el ploteo que nos sirve*/
void SOM::saveWeightsFile(int& i,string file){
	ofstream archi(file.c_str(),ios::trunc);
	double x,y;
	if(archi.is_open()){
		for (unsigned int m=0;m<this->vIcone[i].conec.size();m++){
			x = this->vNeurons[i].weights[0];
			y = this->vNeurons[i].weights[1];
			archi<<x<<"       "<<y<<endl;
			x = this->vNeurons[vIcone[i].conec[m]].weights[0];
			y = this->vNeurons[vIcone[i].conec[m]].weights[1];
			archi<<x<<"       "<<y<<endl;
			//cout<<x<<"  "<<y<<"  "<<y<<endl;
		}
		archi.close();
	}
	else printf("Error de escritura de archivo. \n");
}
void SOM::saveWeightsFile(vector<neurons>& v,string file){
	ofstream archi(file.c_str(),ios::trunc);
	double x,y;
	if(archi.is_open()){
		for (unsigned int m=0;m<v.size();m++){
			x = v[m].weights[0];
			y = v[m].weights[1];
			archi<<x<<"       "<<y<<endl;
		}
		archi.close();
	}
	else printf("Error de escritura de archivo. \n");
}
//-----------------------------------------------------------------------------------------------------------
void SOM::normalize(){
	double xMax = inputs[0].inp[0], yMax = inputs[0].inp[0];
	for (unsigned int i=1;i<this->inputs.size();i++){	
		if (inputs[i].inp[0] > xMax) xMax = inputs[i].inp[0];
		if (inputs[i].inp[0] > yMax) yMax = inputs[i].inp[0];
	}	
}
//-----------------------------------------------------------------------------------------------------------
void SOM::print_patterns(){
	ofstream archi("patterns.txt",ios::trunc);
	if(archi.is_open()){
		for (unsigned int i=0;i<this->inputs.size();i++){
			archi<<inputs[i].inp[0]<<"   "<<inputs[i].inp[1]<<endl;
		}
		archi.close();
		gp("plot [-1:1][-1:1] 'patterns.txt' with points;");
	}
	else printf("Error de escritura de archivo de patrones. \n");
}
//-----------------------------------------------------------------------------------------------------------
void SOM::classification(vector<int>& v){
	int winner,cc1=0;
	double ix,min=99999;
	vClass.resize(this->vNeurons.size());
	for (unsigned int i=0;i<this->inputs.size();i++){
		//Busco que neurona esta más cerca
		for (unsigned int j=0;j<this->vNeurons.size();j++){
			ix=0;
			for (unsigned int k=0;k<this->vNeurons[j].weights.size();k++){
				ix += pow(this->inputs[i].inp[k]-this->vNeurons[j].weights[k],2);
			}
			if(ix<min){
				min = ix;
				winner = j;
			}
		}
		vClass[winner] = v[i]; //clasifico
	}
	for (unsigned int j=0;j<vClass.size();j++) if(vClass[j]) cc1+=1;
	print(1);
	printf("Cantidad de neuronas clase 1: %i\n",cc1);
	printf("Cantidad de neuronas clase 2: %i\n",this->vNeurons.size()-cc1);
	printf("Resumen Neurona con sus respectivas clases. \n");
	for (unsigned int i=0;i<this->vNeurons.size();i++){
		cout<<"Neurona [";
		for (unsigned int j=0;j<this->vNeurons[i].weights.size();j++){
			cout<<this->vNeurons[i].weights[j]<<" ";
		}
		cout<<"] clase ["<<vClass[i]<<"]"<<endl;
	}
}
